function glm = vtc_CreateGLM(hfile, sdm, glmopts)
% VTC::CreateGLM  - create a single study GLM
%
% FORMAT:       glm = vtc.CreateGLM(sdm [, glmopts])
%
% Input fields:
%
%       sdm         SDM object with same number of time points
%       glmopts     optional fields for GLM calculation
%        .threshval threshold value to compute a GLM (default: 100)
%        .threshvol threshold vols (mean value, default: [fix(N/2)+1,N])
%        .tctrans   time course transformation before estimation
%                   'none' or '' no transformation
%                   'z' z-transformation
%                   'psc' %-transformation
%
% Output fields:
%
%       glm         GLM object/structure

% to-do:
%        .arcorr    auto-correlation correction (default: 0)

% Version:  v0.8a
% Build:    9102117
% Date:     Oct-21 2009, 5:36 PM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'vtc') || ...
    numel(sdm) ~= 1 || ...
   ~isBVQXfile(sdm, 'sdm')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end

% get VTC object for some more checks
sc = bvqxfile_getscont(hfile.L);
bc = sc.C;
szvtc = size(bc.VTCData);
numtp = szvtc(1);
szvtc(1) = [];
nlvtc = prod(szvtc);

% further checks on SDM object
scsdm = bvqxfile_getscont(sdm.L);
bcsdm = scsdm.C;
sdmfile = scsdm.F;
if isempty(sdmfile)
    sdmfile = 'Interactive';
end
if bcsdm.NrOfDataPoints ~= numtp
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end

% options
if nargin < 3 || ...
   ~isstruct(glmopts) || ...
    numel(glmopts) ~= 1
    glmopts = struct;
end
if ~isfield(glmopts, 'tctrans') || ...
   ~ischar(glmopts.tctrans) || ...
    numel(glmopts.tctrans) ~= size(glmopts.tctrans, 2) || ...
   ~any(strcmpi(glmopts.tctrans, {'psc', 'z'}))
    glmopts.tctrans = 'none';
else
    glmopts.tctrans = lower(glmopts.tctrans);
end
switch (glmopts.tctrans)
    case {'none'}
        transval = 0;
    case {'psc'}
        transval = 3;
    case {'z'}
        transval = 1;
end
if ~isfield(glmopts, 'threshval') || ...
   ~isa(glmopts.threshval, 'double') || ...
    numel(glmopts.threshval) ~= 1 || ...
    isinf(glmopts.threshval) || ...
    isnan(glmopts.threshval)
    glmopts.threshval = 100;
else
    glmopts.threshval = real(glmopts.threshval);
end
if ~isfield(glmopts, 'threshvol') || ...
   ~isa(glmopts.threshvol, 'double') || ...
    any(isinf(glmopts.threshvol(:)) | isnan(glmopts.threshvol(:)) | ...
     glmopts.thresvol(:) < 1 | glmopts.threshvol(:) > numtp)
    glmopts.threshvol = [fix(numtp / 2) + 1, numtp];
else
    glmopts.threshvol = fix(real(glmopts.threshvol(:)'));
end

% define mask
if ~isempty(glmopts.threshvol)
    tvol = glmopts.threshvol;
    tval = glmopts.threshval;
    maski = false([1, szvtc]);
    for slc = 1:szvtc(3)
        sldat = bc.VTCData(:, :, :, slc);
        maski(1, :, :, slc) = mean(sldat(tvol, :, :), 1) > tval;
    end
else
    maski = true([1, szvtc]);
end

% get iXX
[sdmret{1:2}] = sdm_CalcBetas(sdm, rand(numtp, 1));
iXX = sdmret{2};

% check constant confound
sdmc = bvqxfile_getcont(sdm.L);
psdm = bcsdm.SDMMatrix;
numpred = size(iXX, 2);
if numpred > size(psdm, 2)
    psdm(:, end + 1) = 1;
    sdmc.PredictorNames{end + 1} = 'Constant';
    sdmc.PredictorColors(end + 1, :) = 255;
end

% build GLM structure
glm = BVQXfile('new:glm');
glmc = bvqxfile_getcont(glm.L);

% build complete GLM structure
glmc.ProjectType = 1;
glmc.ProjectTypeRFX = 0;
glmc.FileVersion = 4;
glmc.NrOfSubjects = 1;
glmc.NrOfSubjectPredictors = numpred;
glmc.NrOfTimePoints = size(bc.VTCData, 1);
glmc.NrOfPredictors = numpred;
glmc.NrOfConfounds = 1;
glmc.NrOfStudies = 1;
glmc.NrOfStudiesWithConfounds = 1;
glmc.NrOfConfoundsPerStudy = glmc.NrOfConfounds;
glmc.SeparatePredictors = 0;
glmc.TransformationType = transval;
glmc.Resolution = bc.Resolution;
glmc.SerialCorrelation = 0;
glmc.MeanAR1Pre = 0;
glmc.MeanAR1Post = 0;
glmc.XStart = bc.XStart;
glmc.XEnd = bc.XEnd;
glmc.YStart = bc.YStart;
glmc.YEnd = bc.YEnd;
glmc.ZStart = bc.ZStart;
glmc.ZEnd = bc.ZEnd;
glmc.CortexBasedStatistics = 0;
glmc.NrOfVoxelsForBonfCorrection = sum(maski(:));
glmc.CortexBasedStatisticsMaskFile = '';
glmc.Study = struct;
glmc.Study.NrOfTimePoints = numtp;
glmc.Study.NameOfAnalyzedFile = sc.F;
glmc.Study.NameOfRTCFile = sdmfile;
glmc.Predictor = struct;
for pc = 1:numpred
    glmc.Predictor(pc).Name1 = sprintf('Predictor: %d', pc);
    glmc.Predictor(pc).Name2 = sdmc.PredictorNames{pc};
    glmc.Predictor(pc).RGB = [255, 255, 255; zeros(3, 3)];
end
glmc.Predictor(end).Name1 = 'Constant confound';
glmc.DesignMatrix = psdm;
glmc.GLMData = struct;
glmc.GLMData.MultipleRegressionR = single(zeros(szvtc));
glmc.GLMData.MCorrSS = single(zeros(szvtc));
glmc.GLMData.BetaMaps = single(zeros([szvtc, numpred]));
glmc.GLMData.XY = single(zeros([szvtc, numpred]));
glmc.GLMData.TimeCourseMean = single(zeros(szvtc));

% only work within mask!
maskx = find(maski);

% get masked data
switch transval
    case {0}
        data = double(bc.VTCData(:, maskx));
    case {1}
        data = ztrans(bc.VTCData(:, maskx));
    case {3}
        data = psctrans(bc.VTCData(:, maskx));
end

% calculate betas and PTC
[beta, isdm, ptc] = sdm_CalcBetas(sdm, data);
glmc.iXX = isdm;

% calculate MC R and total SS and set fields
vartc = var(data);
xy = (sdmc.SDMMatrix' * data)';

% create values
glmc.GLMData.MultipleRegressionR(maskx) = std(ptc) ./ sqrt(vartc);
glmc.GLMData.MCorrSS(maskx) = (numtp - 1) .* vartc;
for betac = 1:size(beta, 2)
    glmc.GLMData.BetaMaps((betac - 1) .* nlvtc + maskx) = beta(:, betac);
    glmc.GLMData.XY((betac - 1) .* nlvtc + maskx) = xy(:, betac);
end
glmc.GLMData.TimeCourseMean(maskx) = (1 / numtp) .* sum(data);

% put glm data in global storage
bvqxfile_setcont(glm.L, glmc);
